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NOMENCLATURE 


h Local heat transfer coefficient 

k Roughness parameter 

k f Roughness parameter for natural icing predictions 

k t Roughness parameter for tunnel predictons 

M Mass of water entering control volume 
MVD Median Volume Droplet Diameter, um 

L Latent heat of fusion for liquid water 
LWC Liquid water content, d/m 3 

n Freezing fraction 

S Surface Distance 

T edge Local temperature at edge of boundary layer 

T S urf Local temperature at ice surface 
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EXECUTIVE SUMMARY 


The NASA-Lewis Research Center in Cleveland, Ohio, with the support of the FAA under 
an interagency agreement for icing research, has directed the development of a time- 
dependent, analytical model of the ice accretion process. This model predicts the shape of 
the ice accretion that would form on an arbitrary two-dimensional geometry, such as an 
airfoil, when exposed to icing conditions. The model is embodied in a computer code called 
LEWICE, which has now been operational for several years. Predictions from the code have 
been extensively compared to ice shapes grown under the same conditions in the NASA- 
Lewis Icing Research Tunnel (IRT) and the British RAE tunnel with generally encouraging 
results. In this study, LEWICE ice shape predictions are systematically compared for the 
first time to ice shapes formed on the wing of the NASA-Lewis Twin Otter Icing Research 
Aircraft during research flights in natural icing conditions over several icing seasons. 

The NASA-Lewis Twin Otter is fitted with an array of instrumentation and an on-board tape 
and monitoring system to record the flight parameters and icing cloud properties associated 
with each icing encounter. 

During each flight, a digital tape recorder was used to record the flight parameters, which 
include aircraft speed, altitude, and angle of attack. Icing cloud properties were also 
recorded, including ambient air temperature, liquid water content (LWC), and the droplet size 
distribution, from which the volume median droplet diameter (MVD) is computed. 

An important difference between natural and tunnel ice shape data is the method used to 
document the ice shape profile. Profiles obtained from tunnel experiments are created by 
conventional ice shape tracings. In natural icing, because of the effects due to sublimation 
and erosion, it is generally not feasible to wait until the aircraft returns to the ground to 
obtain a tracing. Therefore, natural ice shape profiles are obtained soon after exiting the 
cloud from a stereo photography system developed jointly by NASA Lewis and the U.S. Air 
Force Arnold Engineering and Development Center (AEDC). This system utilizes two 
Hasselblad 70 mm format cameras mounted in the nose of the Dehavilland DHC-6 Twin 
Otter Icing Research Aircraft. The cameras photograph the ice shape from two different 
viewpoints simultaneously, producing a stereo image which is later photogrammetrically 
analyzed to obtain a profile of the ice shape. 

The main findings of this study are as follows: (1) An equivalent sand grain roughness 
correlation different from that used for LEWICE tunnel comparisons must be employed to 
obtain satisfactory results for flight. (2) Using this correlation and making no other changes 
in the code, the comparisons to ice shapes accreted in flight are in general as good as the 
comparisons to ice shapes accreted in the tunnel. As in the case of tunnel ice shapes, 
agreement is least reliable for large glaze ice shapes at high angles of attack. (3) In some 
cases comparisons can be somewhat improved by utilizing the code so as to take account of 
the variation of parameters such as liquid water content (LWC), which may vary significantly 
in flight. 


v 


PRECEDING PAGE BLANK NOT FILMED 




ANALYTICAL ICE SHAPE PREDICTIONS FOR 
FLIGHT IN NATURAL ICING CONDITIONS 


Brian M. Berkowitz 
Sverdrup Technology, Inc. 

NASA Lewis Research Center Group 
Cleveland, Ohio 44135 

and 

James T. Riley 

Federal Aviation Administration 
Technical Center 

Atlantic City International Airport, New Jersey 08405 


INTRODUCTION 

The NASA-Lewis Research Center in Cleveland, Ohio has directed the development of a 
time-dependent, analytical model of the ice accretion process. This model predicts the shape 
of the ice accretion that would form on an arbitrary two-dimensional geometry, such as an 
airfoil, when exposed to icing conditions. The model is embodied in a computer code called 
LEWICE, which has now been operational for several years. The LEWICE code is part of the 
overall aircraft icing analysis methodology being developed by NASA. Portions of this 
activity, such as the development and evaluation of LEWICE, are being done jointly with the 
FA A. Predictions from the code have been extensively compared to ice shapes grown under 
the same conditions in the NASA-Lewis Icing Research Tunnel (IRT) and the British RAE 
tunnel with generally encouraging results. In this study, LEWICE ice shape predictions are 
systematically compared for the first time to ice shapes formed on the wing of the NASA- 
Lewis Twin Otter Icing Research Aircraft during research flights in natural icing conditions 
over several icing seasons. Such comparisons are essential to the eventual validation of 
LEWICE for specific uses in applied icing analysis by industry and the FAA. 


DESCRIPTION OF CODE 

The NASA-Lewis ice accretion prediction code (LEWICE) is an analytical ice accretion 
model that evaluates the thermodynamics of the freezing process that occurs when 
supercooled droplets impinge on a body (reference 1). LEWICE consists of three main 
modules (figure 1). (1) A flow field calculation employing a Hess and Smith 2D paneling 
method permits the calculation of the air velocity at any point in the flow field. (2) A droplet 
trajectory calculation determines the trajectories of the water droplets driven by the air 
velocities and the resulting impingement pattern of the water droplets on the surface. (3) A 
quasi-steady state surface heat transfer analysis permits calculation of the local ice thickness 
over the geometry, which in turn determines the ice accretion. 

The icing conditions input to the code consist of temperature, LWC, droplet size, velocity, 
and pressure. The surface of the geometry is defined by line segments joining a set of 
discrete body coordinates input by the user. 


1 


The code applies a time-stepping procedure to "grow" the ice, layer by layer. The length of 
each time step is input by the user; they need not all be of the same duration. For the first 
time step, the flow field and droplet impingement pattern are determined for the clean, or 
"uniced," geometry. Next, the heat transfer analysis, leading to a local ice thickness, is 
performed for each segment. The body coordinates are then adjusted to account for the first 
layer of ice accreted. For the second time step, the flow field and droplet impingement 
pattern are recomputed for the "new geometry," i.e., the original geometry with the first 
layer of ice. All the computations of the first time step are now repeated, segment by 
segment, for the iced surface and the body coordinates are adjusted once more to account for 
the addition of the second layer of ice. This procedure is continued until the end of the icing 
growth process has been reached. 

The heat transfer analysis which is applied to each segment is essentially that given by 
Messinger (reference 2), which was a modification of earlier work by Tribus (reference 3). 
A control volume (figure 2) is associated with each segment. The outer boundary of the 
control volume lies just beyond the boundary layer. The lower boundary is the line segment 
itself, which for the first time step forms part of the clean geometry and for later time steps 
forms part of the ice surface accreted up to that time. (Note, therefore, that a control volume 
is always situated on either the clean surface or the ice surface.) A heat balance equation is 
formulated for this control volume which includes terms for all significant processes 
providing heat to or removing heat from the volume. A discussion of all of these terms can 
be found in reference 1 as well as in the paper of MacArthur (reference 4). The dominant 
"heat removal" term, to be discussed shortly, is that for convective cooling; the dominant "heat 
source" term is that for the heat released by the freezing of water entering the control volume. 
The water entering a control volume has two sources: (1) water droplets impinging on the 
surface segment; (2) water "running back" from an adjacent control volume closer to the 
stagnation point. (As will be seen, this "run back" water consists of all water which entered 
the adjacent control volume but did not freeze.) 

Let M denote the mass of all water entering the control volume from either source and let L 
denote the latent heat of fusion for water; then the heat released by freezing of the incoming 
water is 


nML, 

where n, the freezing fraction, denotes the fraction of the incoming water which freezes for 
that segment. If n=l, then all incoming water freezes. If n<l, then a fraction 1-n does not 
freeze. This water will in turn run back into the adjacent control volume further away from 
the stagnation point. 

The discussion thus far shows that, for a given set of icing conditions, each new layer of ice 
is mainly determined by two factors: (1) the impingement pattern calculated for the current 
time step; (2) the segment-by-segment distribution of freezing fractions calculated for the 
current time step, since they determine the runback and the local ice thickness over the 
surface. 

For each segment, the freezing fraction is calculated as follows: Assume that the equilibrium 
temperature, T surf , is 273.15 K. With this assumption, all quantities in the heat balance 
except n can be evaluated. Now solve for n. If the calculation yields a value of n between 
0 and 1 inclusive, the calculation is complete. If n is calculated to be larger than 1, assume 
that n * 1 and that the excess over 1 was because T f is actually smaller than 273.15 K. So 
set n equal to 1 in the heat balance equation, whicn is now solved iteratively (since several 
quantities depend on T surf ) for T f . If n is calculated to be smaller than 0, similar 
reasoning leads to setting n equal to l) and solving iteratively for T surf , which will now be 
larger than 273.15 K. 
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Clearly, a dominant term or terms in the heat balance equation will heavily influence n and, 
through n, the ice shape. As noted earlier, the two dominant terms for nearly all comparison 
cases (which have Mach numbers less than .3) are the term for heat released by freezing of 
incoming water and the term for convective cooling. Since the first includes n but the second 
does not, the second plays a very important role in determining n. The convective cooling 
term has the form 


where h is the local heat transfer coefficient, is the local temperature at the edge of the 
boundary layer, and T surf is the local temperature at the surface. The icing cases most 
challenging to an ice prediction code are the so-called glaze ice cases such that n is 
substantially less than 1 over a significant part of the surface. For any segment such that 0 
< n < 1, T surf = 273.15 by definition. T. is calculated from the local edge velocity Vedge, 
which is provided by the flow field module. Thus calculation of the convective cooling term 
is essentially a matter of calculating h. The method used by LEWICE to calculate h is 
described in the next section. 


HEAT TRANSFER COEFFICIENT AND SURFACE ROUGHNESS PARAMETER 

If the flow field calculation in LEWICE were done using a Navier-Stokes solver, the flow 
field module could provide the heat transfer coefficients h for each surface segment in 
addition to the on- and off-body velocities. Since LEWICE uses a panel method to compute 
an incompressible, potential flow (which requires a much less powerful computer), the values 
for h have to be obtained in another way. LEWICE employs an integral boundary layer 
method for this purpose; the method was described in a 1983 AIAA paper by Mac Arthur 
(reference 4). Makkonen (reference 5) published an interesting study of a similar method in 
1985, discussing the derivation in detail and comparing calculated heat transfer coefficient 
distributions to some of the rather limited experimental data available. Makkonen’s method 
differs from that used in LEWICE in that his method uses a most probable roughness element 
height to compute transition from laminar to turbulent flow in the boundary layer, and an 
equivalent sand grain roughness to compute the heat transfer coefficient in the turbulent 
region. LEWICE uses a single roughness parameter for both calculations. 

In LEWICE, the local value of h depends chiefly on the local velocity at the boundary layer 
edge and the local roughness parameter k. In the icing process, one would expect k to vary 
with both position and time, and ideally LEWICE ought to attempt to account for this. For 
example, for the clean airfoil of the first time step, k might be taken to be constant over the 
entire surface. After each layer of ice is accreted, k might then vary as a function of surface 
position with the nature of the accreted ice (e.g., glaze in the stagnation region, rime further 
aft). It would be desirable if LEWICE could realistically adjust the surface k-distribution as 
the ice accretion develops. Unfortunately, the present state of knowledge is far from making 
this possible. In fact, there seems to be no generally accepted way for quantifying the 
roughness of ice accretions, although some recent work may be promising (reference 6); thus 
no equivalent sand grain roughness values are yet available for iced surfaces. Therefore, at 
this time there is no possibility of treating k in the way described above. 

How, then, is k handled by LEWICE? It uses an empirical correlation for k that was 
developed by Ruff as follows (reference 1). First, a set of icing tunnel ice accretion shapes 
was chosen; for each shape, reliable and complete information on the icing conditions was 
available. The set of ice shapes used were grown in the RAE icing tunnel and are reported 
on by Gent (reference 7). Second, for each experimental shape, a sequence of LEWICE runs 
was done varying the value of k (using a single value for k over the entire surface and for all 
time steps for each run). The value of k which yielded the best agreement between the 
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LEWICE prediction and the experimental shape was identified. At the end of this process, 
a data set had been generated which contained the icing conditions and a "best" value of k for 
each experimental shape. Third, this data set was used to develop a correlation expressing 
k as a function of freestream velocity, LWC, and ambient temperature. 

Since the local h values strongly influence the ice shape, and the local k values in turn 
strongly influence the local h values, it is unfortunate that k cannot be treated in a more 
fundamental way. Unfortunately, as noted above, not enough is known about ice roughness 
to make this possible. However, the following points should be noted. First, many 
comparisons of LEWICE predictions made using this correlation to shapes grown in the IRT 
have been made, with generally encouraging results. The fact that this correlation was 
developed using data from another facility (the RAE tunnel) speaks well for the practical 
utility and applicability of the correlation. Second, once qualitative data on ice roughness 
does become available, it will be useful to relate Ruffs correlation to this data. This could 
be one way of developing a kind of "equivalent ice roughness index." Third, if the same 
correlation gives reasonably good agreement for one facility and generally unsatisfactory 
agreement for another facility or for natural icing in flight, it may point to a fundamental 
physical difference between facilities or between a facility and natural icing in flight. This 
point will also be discussed in more detail later. 


COMMENTS ON THE ICING MODEL 

The icing model described previously relies on a simplified description of the icing process. 
Recent high-speed photo-micrographic films by Olsen (reference 8) have been interpreted 
as showing that fundamental assumptions are in error; of particular importance, it has been 
said that in these films runback is no more than an initial transient phenomenon. Bilanin 
(reference 9) has pointed out that liquid film dynamics, momentum transfer from the air to 
water flow, and other possibly important aspects of the icing process are neglected. Hansman 
and Turnock (reference 10) have proposed a zonal approach to ice roughness based on their 
experimental studies. NASA-Lewis (reference 11) has initiated efforts to "formulate either 
changes to the existing model or an alternate model in order to see if ice shape predictions 
can be improved.” 

Certainly the model can be modified so as to achieve greater physical realism; however, for 
the foreseeable future there will be severe constraints on the accuracy of the predictions that 
any model can achieve. Foremost among these constraints is the accuracy of current icing 
instrumentation. For example, to assume less than a twenty per cent error in LWC or MVD, 
two very important icing parameters, is rather optimistic even in controlled icing facilities 
(reference 12). The performance of LEWICE when compared both to tunnel shapes and 
natural flight shapes must be assessed in light of these constraints, as well as the apparently 
inherently stochastic nature of the aircraft icing process itself. 


RESULTS OF PAST TUNNEL COMPARISONS 

Comparisons of LEWICE predictions to ice shapes grown in the IRT and the British RAE 
tunnel have been encouraging. Although a systematic study has not yet been published for 
tunnel comparisons, NASA studies have shown that LEWICE often provides reasonable 
engineering approximations. 

Examples of comparisons of LEWICE predictions to ice shapes grown in the IRT are shown 
in figure 3. Figure 3a presents a fairly challenging case because of the quantity of ice, the 
type of ice (glaze), and the angle of attack (4 degrees). The ice shape in figure 3b was grown 
at 0 degrees angle of attack; the LEWICE prediction achieves good agreement with the 
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impingement limits and the cross-sectional area, the main discrepancy in shape being in the 
region of the lower horn. These examples are "typical" of tunnel comparisons with LEWICE 
predictions. More "extreme" cases (large amounts of glaze ice, perhaps at high angles of 
attack) can result in poorer comparisons. 


STEREO PHOTOGRAPHY IN NATURAL ICING 

An important difference between natural and tunnel ice shape data is the method used to 
document the ice shape profile. Profiles obtained from tunnel experiments are created by 
conventional ice shape tracings. In natural icing, because of the effects due to sublimation 
and erosion, it is generally not feasible to wait until the aircraft returns to the ground to 
obtain a tracing. Therefore, an alternative method is needed to obtain natural ice shape 
profiles soon after exiting the cloud. 

A stereo photography system developed jointly by NASA Lewis and the U.S. Air Force 
Arnold Engineering and Development Center (AEDC) utilizes two Hasselblad 70 mm format 
cameras mounted in the nose of the Dehavilland DHC-6 Twin Otter Icing Research Aircraft. 
The cameras photograph the ice shape from two different viewpoints simultaneously, 
producing a stereo image which is later photogrammetrically analyzed to obtain a profile of 
the ice shape. Figure 4a shows the layout of the stereo photography system on the Twin 
Otter; detailed descriptions of the stereo photography system and analysis methodology can 
be found in references 13 and 14. One limitation of the system is that the stereo photography 
pairs must be taken in clear air, so that it is not possible to obtain any ice shape data while 
the aircraft is in the clouds. Typically, the ice shape surface can be resolved to within +/- 
0.08 cm, which is thought to be adequate resolution for comparisons between experimental 
and predicted ice shapes. Examples of resultant ice shape profiles produced by this system 
will follow in the results section. 

The interpretation of the stereo photographic ice profiles requires some discussion. The 
profiles for Flights 84-29, 84-34, 85-14, 86-17, and 86-21 are "composite profiles" (reference 
14). That is, rather than representing a cross-section of the ice accretion at a fixed station 
on the wing, the profile includes points within a 4 inch span of the wing (approximately + 
2 inches about a nominal station). Figure 4b illustrates the area on the leading edge that is 
used for stereo photography analysis. A distance of approximately 45 inches is available 
within which a 4 inch section of the wing is used to determine a "composite profile" of the 
ice. Thus an ice profile can be thought of as an overlay or intermingling of several cross- 
sectional ice shapes; a sort of band results and this might be interpreted as a "variability band" 
within which it is hoped that the LEWICE prediction will fall. However, a significant caveat 
must be entered: Some of the points shown may only represent bubbles or cracks below the 
surface of the ice. This would have the effect of moving the inner boundary of the "band" 
closer to the clean airfoil surface. In some cases, there may not be a reasonably well defined 
ice shape; Figure 5 is an example where this is true in the lower region, probably because of 
cracks and bubbles below the ice surface. Such cases were excluded from this study since 
meaningful comparison with prediction would be difficult. 

To better define the ice shape profile, an alternative method of analysis was developed in the 
last year. The profiles produced by this method are cross-sections at a given station on the 
wing. It was possible to obtain profiles of this type for Flights 85-17, 85-24, 86-20, 86-21, 
86-31, 86-33, 87-08, and 87-10. In figure 6, profiles are shown at three different stations 
(at two-inch intervals) for Flight 87-08. These plots are overlaid in figure 7, again yielding 
a kind of "variability band." It would appear that for most attributes the variability is on the 
order of ten per cent. This method generally yields a more well defined ice shape profile. 
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ICING CLOUD PROPERTIES 


The NASA -Lewis Twin Otter Icing Research Aircraft was used for the collection of all 
natural icing data. The aircraft is fitted with an array of instrumentation and an on-board 
tape and monitoring system to record the flight parameters and icing cloud properties 
associated with each icing encounter. 

During each flight, a digital tape recorder was used to record the flight parameters, which 
include aircraft speed, altitude, and angle of attack. Icing cloud properties were also 
recorded, including ambient air temperature, liquid water content (LWC), and the droplet size 
distribution, from which the volume median droplet diameter (MVD) is computed. Three 
liquid water content instruments were operational for most of the flights that were selected 
for this study: 1) Johnson-Williams (J-W), 2) Rosemount, and 3) Leigh. The J-W is a hot 
wire probe while the other two are ice-accretion type instruments (more will be said on this 
shortly). The droplet size distributions were obtained from a Forward Scattering 
Spectrometer Probe (FSSP), a laser instrument that sizes particles and assigns them to one of 
fifteen bins, which ordinarily cover a size range of droplet diameters from 2.0 to 47.0 
microns. For a fuller description of the instrumentation, see reference 15. 

Previous comparisons between instruments in natural icing have shown that the liquid water 
content measurements from the J-W have agreed quite well with the measurements made with 
a rotating multicylinder (RMC) (reference 15). The Rosemount and the Leigh did not agree 
as well, showing much higher data scatter. For this and other reasons, the LWC from the J- 
W was used whenever possible in this study. For several flights the J-W was not available, 
and then either the Rosemount or the Leigh was used. 

A software package developed at NASA-Lewis was used to reduce the data from the digital 
tapes (reference 15). Flight data and cloud properties were averaged at the rate of 1/sec. 
Time-averaged values (averaged over the total icing encounter) of velocity, altitude, angle 
of attack, temperature, LWC and MVD were calculated and are given in Table 1 for each of 
the flights selected for this study. Also given in Table 1 is the total time for the icing 
encounter and the type of ice accreted. For several of the flights, additional processing of 
the flight data was needed and this will be discussed in the results section as warranted. 

Figure 8 presents the complete set of parameter traces for the icing encounter of Flight 86- 
17. Note that although it was a twenty minute encounter, the MVD trace terminates after 
thirteen minutes because of problems with the FSSP. Instrumentation problems of this and 
other kinds were not so rare that such cases could be routinely excluded from the study. In 
this instance, the flight had produced a reasonably "interesting" ice shape for comparison with 
prediction, and was otherwise well-documented. Since the MVD had been very steady for 
the first thirteen minutes, it did not seem unreasonable to assume that it continued around 
the same level for the rest of the flight. Thus the flight was included in the study, using for 
the entire encounter the average MVD computed for the first thirteen minutes. Decisions of 
this kind had to be made for several of the flights included in the study; in a number of cases, 
it was decided that the flight had to be excluded since there was simply too much uncertainty 
concerning the behavior of one or more of the parameters during a large part of the 
encounter. 

Unlike the LWC trace in Figure 8, the trace in Figure 9 is characterized by the presence of 
a significant number of "spikes." Several of the flights had LWC traces of this type. Are 
LEWICE predictions less accurate for flights with "spiky" traces? This question is briefly 
addressed in the results section. 
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MODIFICATION OF THE ROUGHNESS PARAMETER FOR FLIGHT COMPARISONS 


The importance of the roughness parameter k has been stressed. During the first phase of this 
study, Ruff’s tunnel correlation for k (denoted k t ) was used for all predictions. The 
agreement of the predictions with the actual ice shapes was generally not as good as had been 
found earlier for tunnel comparisons. At this point, calculations were made for Flight 86- 
17 using a value of k equal to 25% of k t . There was a marked improvement, as shown in 
figure 10. We then decided to do new runs using a roughness parameter k f = .25k t . In most 
cases, there was improvement, sometimes striking. In no case was there, in our judgement, 
any significant degradation. 

Why the improvement? Consider figure 1 1, which shows the two predictions for Flight 86- 
21, as well as the accompanying plots of h and n for the first time step. Using k f , the 
downward shift in h results in a wider and deeper trough in the plot of n. This causes the 
water to run further back on the first time step before freezing. The situation is similar on 
subsequent time steps; this results in better agreement with experiment for this flight. 

Figure 12 is for Flight 85-14; in this case, there is no significant difference in the predictions. 
The reason is indicated in the plots of h and n for the first time step. Although h again shifts 
downward, the heat balance is such that n remains unity over the entire surface (the change 
in h is reflected in a higher predicted surface temperature in the affected region). Thus the 
predictions for the first time step (not shown) are identical. The freezing fraction remains 
at unity for all four time steps; the slight differences in the final predictions are due to 
secondary effects. 

It must be stressed that the rule k f =.25k t is rather arbitrary. Using a different percentage 
might give somewhat better results. Furthermore, a more sophisticated modification in Ruff’s 
correlation than simply taking a fixed percentage might give improved results. Further 
investigation of this matter is warranted. Nonetheless, this simple rule did result in a definite 
overall improvement in the LEWICE predictions. 

What is the physical explanation for the difference in k for flight predictions versus icing 
tunnel predictions? It could be at least partly due to an actual difference in the nature of the 
ice roughness in the two environments; unfortunately, there is no quantitative data that bears 
on this question. Another hypothesis, for which some independent supporting evidence does 
exist, is that it reflects lower turbulence levels in flight. Hansman and Kirby (reference 16) 
have addressed this matter, and their work indicated that lower turbulence levels in flight 
resulted in reduced heat transfer, as would be expected. Thus the 25% adjustment might be 
regarded as a kind of "turbulence correction." 


RESULTS 

Since the addition of the stereo photography system to the Twin Otter in 1983, stereo 
photography data has been collected over four icing seasons. This data formed the basis for 
this study. An effort was made to establish a set of ice shapes for comparison which, while 
not excluding rime ice accretions, favored glaze and mixed accretions, since these present the 
greatest challenge to an ice prediction code and also entail the worst performance penalties. 
It was also the intent to include flights with a wide range of LWC’s and MVD’s. These 
objectives had to be pursued in the presence of a variety of "practical" considerations. Many 
flights were excluded for one or more of the following reasons. First, in the early years of 
development of the stereo photography system, there were problems of poor image quality 
of the photographs, resulting in difficulties with the analysis of the data; such flights had to 
be excluded from this study. Second, as has been discussed above, in some cases the presence 
of bubbles and cracks below the surface resulted in profiles that could not be used for 
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meaningful comparison with a LEWICE prediction, and such flights were excluded. Third, 
a few of the flights were excluded simply because flight or icing cloud data was not available 
on tape. Fourth, as briefly discussed above, several flights could not be used because of 
instrument failure during the icing encounter. Fifth, some flights were excluded because the 
size of the ice accretion was so small as to render it of little interest for comparison to a 
LEWICE prediction. 

Twelve research flights were finally included in the study. Included in the documentation 
for nearly all the flights was information permitting a classification of the kind of ice; there 
were three rime accretions, five glaze accretions, and four mixed accretions. LEWICE was 
used to predict the ice shapes for each of the twelve flights. The time-averaged values shown 
in Table 1 were used as inputs to the code along with an airfoil chord of 1.98 m. An 
equivalent sand grain roughness (k) was calculated as described earlier based on the time- 
averaged values of velocity, LWC, temperature, MVD, pressure, and angle of attack. In each 
case, four or fewer time steps were used to predict the ice shape. 

Figures 13-24 show the ice shape profiles and the LEWICE ice shape predictions for each of 
the twelve flights. The LEWICE predictions show fairly good agreement with the 
experimental ice shapes except for flights 86-20 and 86-31 (figures 19 and 21, respectively). 
LEWICE tends to predict the impingement limits and cross-sectional area of the ice 
reasonably accurately. The prediction of ice shape is, in general, somewhat less satisfactory. 

The shape of an ice accretion is directly related to the type of ice being accreted. Rime ice, 
which was accreted in Flight 85-14 (figure 15), has a shape which tends to conform to the 
leading edge of the airfoil. For most of the rime ice accretions, LEWICE provides reasonably 
good predictions (as is also true for tunnel comparisons). Glaze ice is often characterized by 
the formation of horns like those shown for Flight 86-20 (figure 19). As noted earlier in the 
discussion of tunnel comparisons, LEWICE sometimes does a relatively poor job of shape 
prediction for "extreme" glaze shapes such as this one; although the shape is poorly predicted, 
the impingement limits and cross-sectional area of the ice are reasonably well predicted. 

All predictions shown in figures 13-24 used as input parameters values that were averaged 
over the entire icing encounter. These parameters all vary to some extent during flight, and 
LWC often varies strikingly. The pattern of variation could theoretically have a dramatic 
effect on the ice shape, as the following purely hypothetical example illustrates. Consider a 
ten-minute icing encounter during which the LWC was relatively steady at 0.5 grams/cubic 
meter and the temperature was always just low enough so that the entire accretion was rime. 
Now consider a second ten-minute encounter during which all parameters were identical to 
the first except that LWC was relatively steady around 0.9 grams/cubic meter during the first 
five minutes and around 0.1 grams/cubic meter during the second five minutes, giving an 
average of 0.5 grams/cubic meter for the encounter. A glaze layer would accrete during the 
first five minutes followed by a rime layer during the second five minutes; the resulting shape 
would be quite different than that for the first encounter, although the average parameter 
values would be the same for both. 

The hypothetical example described in the previous paragraph suggests another potential 
advantage of a time-stepping ice prediction code such as LEWICE. The time steps can simply 
be chosen to coincide with the two five-minute subintervals, and the LWC used for each step 
can be the subinterval average for that step, rather than the overall average for the entire 
encounter. If the heat balance equation is sufficiently accurate, the first predicted layer 
should be glaze and the second rime, so that the final prediction should agree better with 
the actual ice shape than the prediction of a single-time-step code that is forced to use LWC 
averaged over the entire encounter. 

This potential advantage can be an actual advantage only for flights such that the LWC trace 
can be divided into subintervals over which LWC is relatively steady; moreover, the LWC 
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levels associated with the subintervals should differ significantly among themselves. Such 
flights were rare among those we examined. There were only three such that subinterval 
averaging appeared to offer any significant advantage. The prediction using overall average 
LWC is shown in figure 25a, and that using the indicated subinterval averages in figure 25b. 
There appears to be some improvement, but it is not marked; this is a secondary effect when 
compared to the impact of the adjustment in the roughness correlation. 

As noted, several of the flights possessed "spiky" LWC traces. Basing subintervals on the 
spikes is not practical, since far too many time steps would be required; furthermore, there 
is no reason to think it would improve the predictions. Does the use of overall averages for 
these spiky traces tend to result in worse predictions than for the steadier traces? No such 
tendency is revealed when predictions using overall averages for the two classes are compared. 


CONCLUSIONS 

The main findings of this study are as follows: (1) An equivalent sand grain roughness 
correlation different from that used for the tunnel comparisons must be employed to obtain 
satisfactory results for flight. (2) Using this correlation and making no other changes in the 
code, the comparisons to ice shapes accreted in flight are in general as good as the 
comparisons to ice shapes accreted in the tunnel. As in the case of tunnel ice shapes, 
agreement is least reliable for large glaze ice shapes at high angles of attack. (3) In some 
cases comparisons can be somewhat improved by utilizing the code so as to take account of 
the variation of parameters such as liquid water content (LWC), which may vary significantly 
in flight. 


FUTURE EFFORTS 

The NASA icing analysis program can be depicted by the triangle in figure 26, where each 
leg indicates comparison between data obtained by the methods at the adjacent nodes. The 
ultimate goal of the program is the development and validation of computer programs to 
predict ice accretion and attendant performance degradation for flight in natural icing 
conditions. Although this paper is confined to ice accretion shapes, the figure also represents 
concurrent efforts underway for aerodynamic performance codes. The research flights of the 
NASA icing research aircraft are an essential part of the program, and systematic comparison 
of LEWICE predictions (as well as performance code predictions) to flight data will continue 
in the future. Controlled experiments can be carried out in the IRT which are not possible 
in flight, and the comparison of such experiments to code predictions will also continue; this 
activity has played an essential role in LEWICE model development and validation. Finally, 
it is desirable that the IRT simulate icing conditions in flight as realistically as possible, and 
also that researchers understand the unavoidable differences. To this end, a full scale wing 
section for the Twin Otter will be tested in the IRT under conditions as nearly as possible 
duplicating several selected icing flights. 

Efforts already underway to better formulate the fundamental physics underlying the 
LEWICE model (references 8, 9, 10) will continue. To determine if any modifications 
resulting from this research are effective, more accurate measurement of icing conditions, 
particularly droplet size and LWC, is necessary. NASA and the FAA will continue to support 
research in this area as well (Reference 17). 
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Input: 

Body Geometry 
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Method of Solution 


• Flow Field Calculation 

- Hess and Smith 2-D Paneling Method For Potential Flow 

• DropletTrajectory/Impingement Calculation 

- FWG 2-D Particle Trajectory Code 

• Thermodynamic Analysis/lce Accretion Calculation 

- Integral Boundary Layer Equations 
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LEWICE Ice Accretion Prediction Code 
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Illustration of Control Volume Used to Formulate Heat Balance Equation 
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FIGURE 3 

Examples Comparing LEWICE Predictions to Ice Shapes 
Grown in the NASA LEWIS Icing Research Tunnel 
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(b) Reference control points on wing surface and wing fence used in analyzing stereo pair 
photographs. Dimensions are in inches. Control points on wing fence positioned on lines 
that intersect the wing surface at 2% chord spacing. 


FIGURE 4 

Stereo Camera System used to Document Ice Shape 
on Wing Ahead of Wake Survey Probe 
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Figure 6 

Ice Shape Profiles at Three Different Stations 
along Wing Span for Flight 87-08 
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Figure 7 

Overlay of Ice Shape Profiles from Figure 6 For Flight 87-08 





Figure 8 

Parameter Traces for Flight 86-1 7 
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LWC (G/M3) 



Figure 9 

LWC Trace with a Large Number of "Spikes" 
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Figure 11a 

Predicted Ice Shape, HTC, and FF Using Roughness Parameter K t 
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Figure 11b 

Predicted Ice Shape, HTC. and FF Using Roughness Parameter K t 


Figure 11 

Comparison of LEWICE Predictions for Flight 86-21 
Using Roughness Parameters kt ana kf 
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Figure 12 

Comparison of LEWICE Predictions for Flight 85-14 
Using Roughness Parameters kt ana kf 

22 

CffiG*H£L PAGE IS 

OF. POOR QUALITY 









ORIGINAL PAGE f.S 
OF POOR QUALITY 


FIT. 84-29 



■1.0 0.0 1.0 2.0 3.0 


X. INCHES 


Figure 13 

Comparison of LEWICE Prediction to Natural Ice 
Shape for Flight 84-29 
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Figure 14 

Comparison of LEWICE Prediction to Natural Ice 
Shape for Flight 84-34 
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FIT. 85*14 




Figure 15 

Comparison of LEWICE Prediction to Natural Ice 
Shape for Flight 85-14 
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Figure 16 

Comparison of LEWICE Prediction to Natural Ice 
Shape for Flight 85-17 
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FIT. OS-24 
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Figure 17 

Comparison of LEWICE Prediction to Natural Ice 
Shape for Flight 85-24 
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Figure 18 

Comparison of LEWICE Prediction to Natural Ice 
Shape for Flight 86-17 
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FIT. 86-20 
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Figure 19 

Comparison of LEWICE Prediction to Natural Ice 
Shape for Flight 86-20 
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Figure 20 

Comparison of LEWICE Prediction to Natural Ice 
Shape for Flight 86-21 
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FIT. 86*31 




Figure 21 

Comparison of LEWICE Prediction to Natural Ice 
Shape for Flight 86*31 
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Figure 22 

Comparison of LEWICE Prediction to Natural Ice 
Shape for Flight 86-33 
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FIT 87-0* 
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Figure 23 

Comparison of LEWICE Prediction to Natural Ice 
Shape for Flight 87-08 
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Figure 24 

Comparison of LEWICE Prediction to Natural Ice 
Shape for Flight 87-10 
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LEWICE 


Figure 26 

NASA Aircraft Icing Analysis Program 



TABLE 1 - TWIN OTTER TIME-AVERAGED FLIGHT DATA 


FLIGHT 

84-29 

84-34 

85-14 

85-17 

85-24 

86-17 

86-20 

86-21 

86-31 

86-33 

87-08 

87-10 

Duration of 
Icing Encounter, 
min. 

12 

22 

58 

60 

22 

19 

54 

45 

62 

39 


39 

Type of Ice 

Glaze 

Mixed 

Rime 

Rime 

Mixed 

Glaze 

Glaze 

Glaze 

Rime 

Glaze 

Mixed 

Mixed 

Airspeed, m/s 

65.59 

70.37 

71.30 

59.98 

70.88 

69.44 

65.48 

67.70 

70.83 

66.36 

60.80 

66.51 

Static 

Temperature, °C 

- 5.3 

- 6.5 

- 14.7 

- 10.3 

- 14.3 

H 

- 9.5 

H 

- 10.0 

- 5.6 

H 

- 9.2 

Altitude, m 

1648 

1700 

2055 

1267 

2013 

2228 

1878 

1312 

2341 

1206 

1300 

1730 

Liquid Water 
Content, g/m3 
(J.W. unless noted) 

0.25 

0.58 

0.17 

0.23 

0.42 

+ 

0.32 

i 

+ 

0.30 

1 

0.35 

0.26 

0.18 

Median Volume 
Diameter, Microns 

16.0 

14.4 

11.3 

12.4 

15.2 

21.0 

14.9 

14.9 

15.1 

13.3 

13.8 

12.5 

Angle of Attack , 
Degrees 

2.1 

1.6 

1.0 

1 

1.5 

1.6 

1.8 

1.6 

1.2 

1.1 

2.3 

1.7 


+ Leigh 

+ ► Rosemount 
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APPENDIX A 


RATING THE AGREEMENT OF LEWICE PREDICTED ICE SHAPES 
WITH NATURAL ICE SHAPES 

In assessing the results of this study, it is desirable to have criteria of agreement between 
prediction and experiment. The approach suggested here is to consider separately three 
attributes: (1) limits of impingement; (2) cross-sectional area; (3) shape. The first two can be 
treated quantitatively by measuring the percentage difference between prediction and 
experiment, although only a rough visual assessment has been made here. A rating of 1 
roughly corresponds to a difference of 10% or less, 2 to 20%, 3 to 30%, 4 to more than 30%. 
(Note that the 10% difference can be regarded as corresponding to a kind of "variability band" 
reflecting the variability in cross-section along the wing span.) In cases where it was 
especially difficult to distinguish between a rating of 1 and 2, say, a 1 .5 was awarded. The 
third attribute, shape, does not easily lend itself to a qualitative approach. One must ask if 
the prediction captures the important features of the actual ice shape (without predicting any 
additional important features not present). This requires a judgement as to what features 
are important (which may depend on the intended application for the predictions), as well as 
to whether a feature has been "captured" and, if so, how well. We will simply offer our 
admittedly subjective judgement, this time with a rating of 1 corresponding to excellent, 2 
to good, 3 to fair, and 4 to poor, and with halves again allowed in cases of indecision. The 
ratings for all three attributes are given for each flight in table Al. This admittedly crude 
rating system should serve as a stimulus to discussion and perhaps to the formulation of a 
better system; if LEWICE is to be modified to achieve greater physical realism, there should 
be a systematic way to determine if greater predictive accuracy has been achieved as well. 

A summary of the ratings is presented in table A2. Note that the predictions were rated at 
least 2.5 on impingement limits in ten cases out of twelve, on cross-sectional area in ten cases 
out of twelve, and on shape in seven cases out of twelve. 
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TABLE Al. RATINGS OF THE LEWICE PREDICTIONS FOR THE TWELVE FLIGHTS 



Impingement 

Limits 

Cross-sectional 

Area 

Shape 

Flight 84-29 

2 

2 

2 

Flight 84-34 

2.5 

1 

2 

Flight 85-14 

1 

1.5 

1.5 

Flight 85-17 

3 

1.5 

3 

Flight 85-24 

1.5 

1.5 

2.5 

Flight 86-17 

1 

1.5 

1 

Flight 86-20 

1.5 

1.5 

4 

Flight 86-21 

1.5 

1.5 

2.5 

Flight 86-31 

3.5 

4 

3 

Flight 86-33 

2 

2 

3.5 

Flight 87-08 

1.5 

2 

2.5 

Flight 87-10 

2.5 

3 

3 

Rating 

TABLE A2. SUMMARY OF RATINGS 

Impingement Cross-sectional 

Limits Area 

Shape 

1 

2 

1 

1 

1.5 

4 

6 

1 

2 

2 

3 

2 

2.5 

2 

0 

3 

3 

1 

1 

3 

3.5 

1 

0 

1 

4 

0 

1 

1 
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